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Abstract 

We report results for the interaction measure, pressure and energy density for nonzero temperature QCD 
with 2+1 flavors of improved staggered quarks. In our simulations we use a Symanzik improved gauge 
action and the Asqtad 0{a^) improved staggered quark action for lattices with temporal extent Nt =4 and 6. 
The heavy quark mass is fixed at approximately the physical strange quark mass and the two degenerate 
Ught quarks have masses ruud ~ 0. 1 or 0.2ms. The calculation of the thermodynamic observables employs 
the integral method where energy density and pressure are obtained by integration over the interaction 
measure. 

PACS numbers: 12.38 Gc, 12.38 Mh, 25.75 Nq 
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I. INTRODUCTION 



Ordinary hadronic matter undergoes a qualitative change into a quark-gluon plasma (QGP) at 
high temperatures and/or densities. The QGP is a new state of strongly interacting matter in which 
the basic constituents, quarks and gluons, are "freed" from the color confinement of low tempera- 
ture hadrons. The phenomenon of color confinement is attributed to the non-perturbative structure 
of the QCD vacuum at zero temperature. At high temperatures (and/or densities) this picture is 
modified to allow a deconfining transition. However, the character of the QGP at temperatures up 
to at least several times the transition temperature {Tc ~ 170 MeV) remains non-perturbative, since 
in this temperature range the strong coupling constant is still of 0(1), and the fundamental degrees 
of freedom are more complex than simply free quarks and gluons. Currently lattice QCD is the 
only theoretical tool that is suitable for tackling this inherently strongly coupled system from first 
principles. 

The QGP is studied experimentally in heavy-ion collisions at RHIC and CERN, in which the 
accessible temperature range is up to about 3Tc yj]. The data from these experiments are mostly 
interpreted through hydrodynamical models [|2|], which take the equation of state (EOS) of the low- 
and high-temperature phases as essential inputs. The hydrodynamic models that include the QGP 
as the high-temperature phase use an ideal gas EOS for quarks and gluons which, considering 
the temperature range, is bound to be an unsatisfactory approximation, perhaps accounting in 
part for discrepancies between some of the current predictions and the experimental data. This 
difficulty can be addressed by a realistic lattice QCD calculation of the EOS to serve as input for 
the hydrodynamics equations. 

The importance of a realistic EOS of the QGP is not limited to the heavy-ion experiments. The 
EOS is also relevant to cosmology, since it is believed that the QGP existed microseconds after 
the big bang. For example, the relic density of weakly interacting massive particles is sensitive 
to the EOS of the QGP at these early stages of the formation of the Universe yD. Another area 
of potential application of the EOS is in the study of phenomena in the interior of dense neutron 
stars, where again the QGP is likely to exist. 

The determination of the EOS through numerical simulation of lattice QCD is challenging, 
since it requires a precise determination of differences between high and low temperature quan- 
tities that have inherent ultraviolet divergences. Thus the most extensive simulations to date are 
carried out on rather coarse lattices (Nr = 4) Improving the gauge and fermion actions 
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1^ 7|] helps reduce lattice artifacts as does decreasing the lattice spacingjfM = 6) [|7Ll8|]. It is also 
important to carry out simulations with a realistic light quark spectrum ^[M. 

In this paper, we report results of a simulation of the QCD EOS at Nf = 6 with 2 + 1 light flavors 
of 0{a^) tadpole-improved (Asqtad) staggered quarks. The gauge action we use is a Symanzik 
0{cP') tadpole-improved one as well. Preliminary accounts were given at the Lattice 2005 and 
2006 conferences [7, ^. The inclusion of the strange quark is of interest to the phenomenological 
studies of the QGP since it can change the order of the phase transition and influences strangeness 
production in the heavy-ion experiments. To determine the EOS we use the integral method where 
the pressure and the energy density are calculated through an integration over the interaction mea- 
sure [10]. The paths of the integration in the bare parameter space are approximately trajectories 
of constant physics. Along a trajectory of constant physics the heavy quark mass (m^) would be 
fixed to the physical strange quark mass and the m^/mp ratio would be kept constant. We ap- 
proximate two such trajectories for which m,i/mp ^ 0.3 and 0.4, which correspond to light quark 
masses niud ~ O.lm^. and O.lnis, respectively. Our calculations are performed dXNt = 6 for both 
trajectories, and we have an additional Nt = 4 result for the m^^ ^ 0. 1 trajectory. In this work 
we compare the EOS obtained using, first, the data from the two different trajectories and, second, 
from the data with different Nf. We find that the differences are small in both cases. 

II. THE INTEGRAL METHOD FOR THE EOS DETERMINATION 

In this section we give a brief description of the formalism of the integral method as applied 
to the specific improved actions that we use. The analytic form of the EOS is derived from the 
following thermodynamics identities 

p dlnZ 



ainz 

eV = 



d{l/T) 



T dV 



InZ , ^ T dlnZ 
J V V dma 



where e is the energy density, p is the pressure and / is the interaction measure. The spatial 
volume is y = N^a^ for lattice spacing a, and the temperature is T = 1 / (Nta) . The derivative 
of the partition function Z with respect to the logarithm of the lattice spacing. In a, should be 
understood as taken along a trajectory of constant physics. In the explicit form of the partition 
function 

Z = y^J[/exp|-5g + £(/i//4)Trln[M(am/,t/,ao)]| (2) 
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the gauge action is given by Sg = Spi + Sn + Spg, with 

^pi = p E (i-z'^v) 

X,IJ<V 

Sn = Prt £ (l-^,.v) 

X,fJ<V 
X,/J<V<CT 

The real parts of the traces of the 1x1 plaquette - P^y, the 1x2 and 2x1 rectangle sum - R/jy, and 
the 1 X 1 X 1 parallelogram - C^vo> are all divided by the number of colors. The gauge couplings 
in the above are defined as 

P = 10// 

Prt = --^{l+0AS05a,) (3) 
20mq 

P = -4o. 03325a, 

fora5 = — 41n(Mo)/3.0684 and uo= (P)^^^- The fermion matrix M(am/,t/, mo) corresponds to the 
Asqtad staggered action for a specific flavor /. 

Using the identities in Eq. ([T]) and the explicit form of Z we obtain the EOS expressions 

7^4 ^ _^JLa(^p) _ U^A{R) - 16^A(C) (4) 
Jlna dlna dlna 



f ^ 



d(mfa) , ,_ , duQ , /_ dM 
dlna dma \ duQ 



pa^ = - r" I{a'){aYdlna' (5) 

Jlnao 

za^ = (/ + 3p)a^ (6) 

The various fermionic and gluonic observables in the EOS are calculated at nonzero temperature 
(fixed Nt < Ns) and on zero-temperature lattices (Nt > Ns). The symbol A in the EOS expressions 
denotes their nonzero and zero-temperature differences. All measurements are taken along a tra- 
jectory of constant physics, which we parameterize with the lattice spacing a. The couplings p, 
Prt, Ppg, masses anif and tadpole factor uq are all functions of In a along this trajectory. We use 
these functions to determine the derivatives of the bare parameters with respect to Ina as needed 
for the EOS. The lower integration limit, Inao, in Eq. ([5]) should be taken at a coarse enough lattice 
spacing that the pressure difference is negligible. 
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III. SIMULATIONS ON TRAJECTORIES OF CONSTANT PHYSICS 



As already mentioned in the previous sections, for our simulations we use the one-loop 



IJ. Both 



Symanzik improved (Liischer-Weisz) gauge action and the Asqtad quark action [|lll. ll2l 
actions are tadpole improved and the leading discretization errors are 0{a^<Xs, • There are many 
features of the Asqtad action that make it well suited for high temperature studies. It has excellent 
scaling properties leading to faster convergence to the continuum limit and the dispersion relations 
for free quarks are much better than the ones for the standard Wilson or staggered actions. An- 
other very important property of the Asqtad action is the much reduced taste symmetry breaking 
compared with the conventional staggered action. All this translates into decreased lattice artifacts 
above the phase transition. 

It would seem important for studying the strange quark physics of the plasma that the kaon 
mass be heavier than the pion. In the staggered fermion scheme each meson state appears in 
a taste multiplet of 16. With improvement of the fermion action the splittings are considerably 
reduced. The splittings in meson mass squared are expected to vanish in the continuum limit as 
a^Ctjy. Shown in a log-log plot in Fig.[T]are pion taste splittings relative to the Goldstone pion mass 
for five lattice spacings. The solid line shows the expected scaling slope. The trend is consistent 
with the scaling prediction. Shown also are splittings of the lowest member of the kaon multiplet, 
relative to the Goldstone pion mass for the two choices of ntud/ms in the thermodynamics study. 
The vertical lines locate the lattice spacing at the crossover temperature Tc (about 190 MeV for 
our unphysical light quark masses) for various Nt. Note that the temperature then increases as we 
move to the left. Our nonzero temperature studies are at A^;^ = 4 and 6. hi Nt = A the lightest 
kaon at has approximately the same mass as the lowest non-Goldstone pion. As the figure 
shows, dX Nt = A the kaon and pion taste multiplets are non-overlapping at approximately T > 2Tc. 
At the Nt = 6 crossover the situation has improved and the multiplets are non-overlapping at 
approximately T > 4Tc/3. Clearly, Nt = S would be even better for this action. 

AtNt = 6 the taste- splitting is about half as large as at N = 4. One of our goals was to determine 
to what extent the increase in N from 4 to 6 influences the EOS. 

In our simulations we use the dynamical R algorithm [llSll with step size equal to the minimum 
of 0.02 and 2am„^/3. For some runs the step size was chosen to be even smaller. Our aim is 
to generate zero- and nonzero-temperature ensembles of lattices with action parameters chosen 
so that a trajectory of constant physics (m,i/mp = const) is approximated. Along the trajectory 
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FIG. 1: Pion taste splitting relative to the Goldstone pion mass in units of ri = 0.318(7) (4) fm vs. the 
lattice scaling variable (a/ri)^av(a)^ in a log-log plot. Here we take av{a) = 127l/[541n[(3.33/<3A)] with 
A = 319 MeV. The rising line has slope 1. The fancy diamonds locate the kaon splittings {mj^ — my)r\ for 
^ud ~ 0.2ms. The fancy crosses do the same for niud ~ 0.1 m^. The vertical lines indicate the approximate 



lattice spacing at the crossover temperature for various Nf. Data are from lll4ll and unpublished simulation 



results. The pion taste assignments are given in the gamma matrix basis. The taste singlet is denoted 71^. 

the heavy quark mass is tuned to the strange quark mass within 20%. We work with two such 
trajectories: niud ~ O.lnis (m,i/mp ^ 0.4) and m^d ~ 0. 1 (m^i/mp ^ 0.3) as shown in Fig. |2l 

The construction of each trajectory begins with "anchor points" in (3, where the hadron spec- 
trum has been previously studied and the lattice strange quark mass has been tuned to approximate 
the correct strange hadron spectrum lll4ll . We adjusted the value of aniud at the anchor points to 
give a constant (unphysical) ratio m^t/mp. Between these points the trajectory is then interpolated, 
using a one-loop renormalization-group-inspired formula. That is, we interpolate ln(am5) and 
\n{amud) linearly in p. Since we have three anchor points for the rriud ~ O.lnis trajectory, namely 
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FIG. 2: Plot of the two trajectories of constant physics in the (am„^,P) plane. 

P = 6.467, 6.76, and 7.092, our interpolation is piecewise linear. For the trajectory niud ~ 0. 1 we 
use two anchor points at (3 = 6.458 and 6.76. Explicitly the parameterization of the niud ~ 0.2ms 
trajectory is 

, ln(0.050/0.0820) \ 



amv 



0.082exp (^(P- 6.4674) j, Pg [6.467,6.76] 
0.05 exp((P- 6.76) i^gg^), p G [6.76,7.092] 

0.01675exp((P-6.4674)if«^), Pg [6.467,6.76] 
0.01 exp((P- 6.76) iffg^), P G [6.76,7.092]. 

The parameterization of the mud ~ 0. 1 m^ trajectory with P G [6.458, 6.76] is 



(7) 



(8) 



am, 



^/o . ,Jn(0.082/0.05) 

ln(0.0082/0.005 



am^d = 0.005 exp (P-6.76) 



(9) 
(10) 



(6.458-6.76) 

For both trajectories, for values of P out of the interpolation intervals, the parameterization formu- 
lae are used to perform extrapolations. The run parameters of the two trajectories at different Nt 
are summarized in Tables HI Ull and Hill It is apparent from Eqs. (7) and (8) that values of the quark 
mass ratio m^d/ms in Table HlI] deviate slightly from 0.2, since it was our initial intention to keep 
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the hadron masses constant instead. In subsequent practice, as in the O-lm^ trajectory, we chose 
the more convenient alternative of keeping the quark mass ratio fixed. 

For the purpose of the EOS determination, the trajectories of constant physics are most conve- 
niently parameterized by the lattice spacing a, as discussed at the end of the previous section. To 
calculate the various derivatives of the bare parameters with respect to In a we need to determine 



the functional dependence lna((3) . The lattice spacing is determined using the method of Ref. [|l4ll . 
On a large set of zero-temperature ensembles the static potential is measured to determine the mod- 
ified Sommer parameter ri jl^ in lattice units. Specifically, ri is defined by r^F^q{ri) = 1. All 
able measurements of ri/a are then fit to the following asymptotic-freedom-inspired form 
11 



aval 



I117L 



a _ co/(g')+C2rr(g')+C4gV(g") 



(11) 



The definition of 



f(g^) = ^tog^yb,/i2tl)^-m2bos') (12) 



involves the universal beta-function coefficients for massless three-flavor QCD, bo and bi. The 
coefficients cq, C2 and C4 are 

CO = coo + {coiuamud + coisams)/f{g^)+co2{2amud + ams)^/f^{g^) 
C2 = C2o + C2i(2ami,^ + am,.)//(^^) 

C4 = C40 
d2 = d20, 

where cqq = 46.766(447), cqu, = 0.526(122), coi, = 0.1817(708), C02 = -0.00403(204), C20 = 
-4.702(175) X 10^, C21 = 3.321(511) x 10^, C40 = 3.943(84) x 10^ and J20 = 1.276(484) x 10^. 
The fit has %^/DoF ^1.3 and a confidence level of approximately 0.13. This parameterization 
provides a determination of r 1 /a along our trajectories of constant physics (Fig. [3]) . Independently 
of the fit, the absolute scale for a is set from a determination of the T(25 —IS) splitting on a subset 
of these zero-temperature ensembles 2^. An extrapolation to zero lattice spacing then gives 



ri = 0.318 (7) (4) fm [14]. This value was used in conjunction with the above parameterized value 
of ri/a to define the physical lattice spacing in our simulations. 
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0.012 ' — ' — ' — ' — ' — ^ — ' — ' — ' — ' — ^ — ' — ' — ' — ' — ^ — ' 
6.0 6.5 7.0 7.5 

FIG. 3: Inverse lattice spacing in units of f{g^)ri vs. gauge coupling (3 = based on the best fit 

parameterization Eq. ([TT]) . The fitting function is evaluated along the two lines of constant physics, namely 
JT^ud ~ 0. 1ms and 0.2m^ . It is derived from forty measured values of ri/a. Eight of them lie on the trajectories 
of constant physics and are plotted here. 

IV. EQUATION OF STATE RESULTS AND CONCLUSIONS 

In the previous sections, we have outlined the method we follow to determine the temperature 
dependence of the bulk thermodynamic quantities, namely, the interaction measure, pressure and 
energy density, which constitute the EOS for the quark-gluon system. In this section we present 
our numerical results. 

According to the integral method, at the base of our calculation is the determination of the 
interaction measure, which is straightforward from Eq. dH). The nonzero temperature value of 
the interaction measure needs to be corrected for the zero-temperature contributions. This cor- 
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FIG. 4: The interaction measure is shown for both of the trajectories of constant physics and the different 
N,. 

rection is done for about half of the runs by directly measuring the zero-temperature values of 
the fermionic and gluonic observables involved in Eq. dH) and subtracting their resultant zero- 
temperature contribution from the interaction measure at nonzero temperature. For the rest of the 
runs, the zero-temperature correction is calculated by making local interpolations. We need to 
determine as well the derivatives d^/d\na, d^n/dlna, d^^g/dlna, d{mfa)/d\na and duo/d\na 
for each trajectory. For this purpose we take derivatives of the ln(amj,^) and \n{ams) trajectory 
parameterizations, using Eq. ([7]) - Eq. (flOl) . polynomial fits to mo(P) for both trajectories, and the 
a/ri fit from Eq. (fTTI) . Figure |4] shows the interaction measure as a function of the temperature for 
both trajectories of constant physics and A^/s. 

The pressure is obtained from the interaction measure by integration (Eq. ([5])) using the trape- 
zoid method and Fig. \5\ shows our results. Using both the interaction measure and the pressure, 
we calculate the energy density (Eq. Q). The results are presented in Fig.[6l The statistical errors 
on all of the thermodynamic quantities are calculated using the jackknife method, and we ignore 
the insignificant errors on the derivatives of the bare parameters with respect to the lattice scale 
mentioned above. The EOS data in the Figs.|4]-[6]is corrected for the systematic errors due to the 
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FIG. 5: The temperature dependence of the pressure for both of the trajectories of constant physics and the 
different Nf. The continuum Stefan-Boltzmann limit for 3 massless flavors is also shown. 

finite step size, which are discussed later in this section and in more detail in the Appendix, and 
the choice of the lower integration limit in Eq. ([5]). 

From our results for the EOS we find that at the highest studied temperature (~ 380 MeV for 
Nt = 6 and ~ 570 MeV for Nt = 4 ) the energy density is about 10 - 15% below the Stefan- 
Boltzmann three-flavor limit, which is evidence that strong interactions between the plasma con- 
stituents persist in the high temperature phase at several times Tc. The comparison of the EOS for 
the two trajectories of constant physics at Nt = 6 shows some small differences. There is a differ- 
ence in the interaction measure maxima, with the one from the ^ 0.2 ms trajectory somewhat 
larger than the O.lm^ trajectory one. Also, the pressure on the niud ^ 0.2 ms trajectory at coarse 
lattice spacing (see Fig. |5]) is gradually becoming slightly larger with temperature than that on the 
0. 1 ms trajectory. This result is contrary to expectation. We consider that the cause is the accumu- 
lation of various systematic errors in the pressure calculation which we discuss later in this section. 
As a whole, the reduction of the mass of the degenerate light quarks from m^d ~ 0.2 to 0. 1 
does not affect dramatically the basic thermodynamic properties of the system. We see a lot of 
similarity between the EOS at A^^ = 4 and 6 for the mud ~ 0. 1 trajectory. The main differences 
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FIG. 6: The temperature dependence of the energy density for both of the trajectories of constant physics 
and the different Nf. To faciUtate comparison with the ideal gas case, the continuum Stefan-Boltzmann Umit 
for 3 massless flavors is also shown. 

between the two available Nt results is again in the interaction measure, where the maximum at 
A^f = 6 is higher. Although the discretization artifacts at A'^ = 4 are known to be larger than in the 
Nt = 6 case, we find that their effect on the EOS is not very pronounced. 

Our EOS calculation can be affected by the following systematic errors: finite volume effects, 
finite step-size effects, the error in the determination of the lower integration limit in Eq. ([5]), 
possible deviations from the trajectories of constant physics and the uncertainties in the scale 
determination from Eq. (fTT)) . First we discuss the scale determination error. The lattice spacing 
is most difficult to obtain for the Nt = 4 case in the low temperature region, where the lattices are 
coarse. We have estimated that a 5% error on the lattice scale there gives up to a three-sigma effect 
in the region of the interaction measure maximum. This translates as well into up to two-sigma 
effects on the energy density and pressure at high temperatures, since errors accumulate in the 
integration needed to obtain these quantities. 

To address the question of the finite volume effects we have conducted a set of runs at N = 4 
with parameters from the niud ~ 0.1 trajectory on lattices with smaller spatial volume - Vj = 8^. 
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FIG. 7: Volume dependence of results for Nt = 4 with ruud ~ O.lm^. Empty symbols are used for small 
volume (Vs = S^) results, and filled symbols are used for large volume (V^ = 12^ or 16^) results. The energy 
density, pressure and interaction measure are plotted using diamonds, squares and circles, respectively. The 
data are not corrected for any systematic errors. We see no statistically significant volume dependence. 

In Fig.|7]the thermodynamic quantities calculated using = 8^ are compared with the ones on the 
larger spatial volumes - Vy = 12^, 16^. We find no statistically significant difference which leads 
us to conclude that in our calculation the finite volume effects are negligible. 

The determination of the lower integration limit in Eq. dS]) is potentially another source of 
systematic error. The lowest available temperature in our calculations is around 135 MeV at 
Nt = 4, and 149 MeV at A^^ = 6. To estimate the pressure at these temperatures we calculate the 
pressure of an ideal Bose gas of pions with masses similar to those in our simulations. The true 
Goldstone pions on the physics trajectories have mass of ~ 270 MeV and the rest of the members 
of the pion multiplet are heavier. We estimated the heavy pion masses using extrapolations of 
available data for the taste splitting in the pion multiplet, summarized in Fig. [T] Including all of 
the pions according to their degeneracy, we estimate p/T^{T = 135 MeV) ~ 0.02 and p/T^{T ~ 
149 MeV) ~ 0.03 with about 30% uncertainty in these values. Both of these estimations are 
comparable or a bit larger than the size of the statistical error on the pressure at the closest available 
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low temperatures. Consequently we have corrected the pressure and energy density by adding 
them to the data. At high temperatures this correction is smaller than the statistical error. 

The error due to deviations from the trajectories of constant physics would be largest in the 
Nt = A case, for which the points around the transition region and at lower temperature were 
obtained by extrapolations using Eqs. (9) and (10). Indeed, a later spectrum calculation near the 
transition, at (3 = 6.2, showed that there is about a 10% difference from the target value for mji/mp. 
However, considering that the differences between the two Nt = 6 trajectories, for which mji/mp 
differs by about 30%, is no more than four sigma, we estimate the effect at about one sigma in the 
transition region and smaller outside of it. 




f i « * ? , \ , \ , \ , 

150 200 250 300 350 400 

T(MeV) 



FIG. 8: Effect of step-size corrections for Nt = 6 with O-lniy. We use filled (open) symbols to plot 

uncorrected (corrected) results. The symbols for the energy density, pressure and interaction measure are 
diamonds, squares and circles, respectively. 

The last potentially significant source of systematic error is the finite step size used in the R 
algorithm. For Nt = A and 6 we have carried out a set of test simulations at a larger step size in the 
R algorithm to estimate their effect. In addition we have performed some RHMC 3] calculations 
to complement our finite step-size study. Our analysis of the results is presented in the Appendix. 
We find that the effect of the step-size corrections to the gauge observables on the interaction 
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measure is no larger than the size of our statistical error along the m„j ^ 0. 1 trajectory and 
negligible along the mud ~ 0.2 ms one. The effects of the correction to the fermionic observables 
for both trajectories is small enough to be ignored. We use the empirical formula in Eq. (lAll) with 
the parameters in Eq. (IA2I) to compute the correction to the three gauge loop observables for the 
i^ud ~ O.lm, trajectory only. We do not correct the m^d ~ O.lm^ trajectory for finite step-size 
effects due to their smallness. Figure [8] shows the EOS for the mud ~ 0. 1 m.,, Nt = 6 case with the 
finite step-size correction compared with the uncorrected case. The correction is no larger than 
our statistical errors. As explained in the appendix, we estimate the uncertainty in the correction 
itself to be about 50%. 

ACKNOWLEDGMENTS 

We thank Frithjof Karsch and Peter Petreczky for useful discussions on the scale determination. 
Computations for this work were performed at Florida State University, Fermi National Accelera- 
tor Laboratory (FNAL), Indiana University, the National Center for Supercomputer Applications 
(NCSA), the National Energy Resources Supercomputer Center (NERSC), the University of Utah 
(CHPC) and the University of California, Santa Barbara (CNSI). This work was supported by 
the U.S. Department of Energy under contracts DE-FG02-91ER-40628, DE-FG02-91ER-40661 
and DE-FG03-95ER-40906, and National Science Foundation grants PHY05-55243, and PHY04- 
56556. 



[1] K. Adcox et al., Nucl.Phys. A 757, 184 (2005), [nucl-ex/04 10003]. 

[2] P. Huovinen and R V. Ruuskanen (2006), [nucl-th/0605008]. 

[3] M. Hindmarsh and O. Philipsen, Phys. Rev. D71, 087302 (2005), [hep-ph/0501232]. 

[4] T. Blum, L. Karkkainen, D. Toussaint, and S. A. Gottlieb, Phys. Rev. D51, 5153 (1995), [hep- 
lat/94 10014]. 

[5] F. Karsch, E. Laermann, and A. Peikert, Phys. Lett. B478, 447 (2000), [hep-lat/0002003]. 

[6] J. Engels et al., Phys. Lett. B396, 210 (1997), [hep-lat/9612018]. 

[7] C. Bernard et al., PoS LAT2005, 156 (2005), [hep-lat/0509053]. 

[8] Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo, JHEP 01, 089 (2006), [hep-lat/05 10084]. 



16 



[9] C. Bernard et al., PoS LAT2006, 139 (2006), [hep-lat/0610017]. 
[10] J. Engels, J. Fingberg, F. Karsch, D. Miller, and M. Weber, Phys. Lett. B252, 625 (1990). 
[11] K. Orginos and D. Toussaint (MILC), Phys. Rev. D59, 014501 (1998), [hep-lat/9805009]. 
[12] D. Toussaint and K. Orginos (MILC), Nucl. Phys. Proc. Suppl. 73, 909 (1999), [hep-lat/9809148]. 
[13] G. P Lepage, Phys. Rev. D59, 074502 (1999), [hep-lat/9809157]. 
[14] C. Aubin et al., Phys. Rev. D70, 094505 (2004), [hep-lat/0402030]. 

[15] S. A. Gotdieb, W. Liu, D. Toussaint, R. L. Renken, and R. L. Sugar, Phys. Rev. D35, 2531 (1987). 
[16] C. W. Bernard et al., Phys. Rev. D62, 034503 (2000), [hep-lat/0002028]. 
[17] C. R. Allton (1996), [hep-lat/9610016]. 

[18] C. R. Allton, Nucl. Phys. Proc. Suppl. 53, 867 (1997), [hep-lat/9610014]. 

[19] M. Wingate, C. T. H. Davies, A. Gray, G. P Lepage, and J. Shigemitsu, Phys. Rev. Lett. 92, 162001 

(2004), [hep-ph/0311130]. 
[20] A. Gray et al., Phys. Rev. D72, 094507 (2005), [hep-lat/0507013]. 

[21] M. A. Clark, A. D. Kennedy, and Z. Sroczynski, Nucl. Phys. Proc. Suppl. 140, 835 (2005), [hep- 
lat/0409133]. 

[22] C. W. Bernard et al. (MILC), Phys. Rev. D55, 6861 (1997), [hep-lat/96 12025]. 



APPENDIX A: STEP-SIZE DEPENDENCE OF EOS OBSERVABLES 



With the R algorithm the integration step size in the molecular dynamics evolution must be 
chosen small enough to achieve the desired accuracy in observables of interest. For most practical 
purposes we have found errors at our standard small production step sizes to be insignificant. (For 
a recent test see Table IX of lll4ll .) However, the observables required for the equation of state 
must be measured to a very high accuracy, since the small differences between the hot and cold 
measurements are sensitive to even small systematic errors. The most important observable in this 
regard is the plaquette. For the present study we have developed a rough empirical method for 
estimating and correcting for these errors in our simulations. 

To estimate the step- size error within the R algorithm requires carrying out simulations at a 
range of step sizes and determining the change in the observable as the step size tends to zero. 
We have carried out a number of such tests on hot and cold Asqtad lattices, measuring most of 
the observables needed for the equation of state. The RHJVIC algorithm, which we incorporated 
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into our code at the end of this study, does not suffer from such step size errors. Thus for the 
purpose of modeling step-size corrections we include results of some RHMC calculations. For 
our previous study of the equation of state with the unimproved gauge action and naive staggered 
fermion action, we made extensive measurements of the step-size dependence of the plaquette and 
chiral condensate [|22n . 

The leapfrog-inspired R algorithm is specifically designed to be a second order integration al- 
gorithm. That is, the truncation error at the end of a molecular dynamics trajectory of fixed length 
decreases with the integration step size £ as In Figs. [9] and [10] we show the step-size depen- 
dence of the plaquette and chiral condensate for the improved action for one pair of ensembles. 
On the larger-volume zero-temperature lattices, we have found that the variation of the plaquette 
with decreasing step size shows more apparent curvature over this range of step sizes than do the 
smaller volume high temperature lattices. Consequently, as shown, we fit the low temperature 
results to a quadratic in Clearly both low and high temperature values are subject to correction. 
The corrections tend to cancel in the difference. For the improved action the slopes for all seven 
observables needed for the equation of state are tabulated in Tables HVl and IVl For the high temper- 
ature ensembles the slope is determined by a linear fit in For the zero temperature ensembles 
it is determined from results at our smallest available pair of values of £^, treating an RHMC step 
size as 0, of course. 



1.606 



□ 1.602 



1.600 



= 0. Ims 




0.09 



0.08 



0.07 



0.06 



0.05 I — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' — ' 
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m^^=0. 1ms 



X 16 x48 
o 16^x6 



FIG. 9: Plaquette (left panel) and chiral condensate (right panel) vs the squared step size £^ for the improved 
action for the ensemble at p = 6.458, am„^ = 0.0082 and anis = 0.082. The squared step size for the 
production of this ensemble is 0.00003 



The step-size correction depends largely on the size of the fermion force, which is computed 
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FIG. 10: Same as Fig.|9]but for the ensemble at P = 6.467, amud = 0.01676 and ants = 0.0821. The squared 
step size for the production of this ensemble is 0.0001 

by inverting the fermion matrix. Small eigenvalues dominate the inverse. The smallest eigenvalue 
is controlled by the light quark mass. Since the chiral condensate is also determined from the 
inverse of the fermion matrix, we would expect the chiral condensate and light quark mass to be 
natural parameters for the step size error regardless of temperature. Consequently, along a chosen 
trajectory of constant physics, we parameterize the step-size slope of the plaquette for both high 
and low temperature ensembles as a polynomial in the chiral condensate. For the m^^ ^ 0. 1 
trajectory we find it is modeled reasonably well by the following quadratic form as shown in 
Fig.im 



dP 



b{x-OAf + m{x-OA)+c, 



(Al) 



where x = {^^)ud is the light quark chiral condensate in lattice units. 
The best fit values are 

b = 250(177), m = -109(54), c = -6.1(2.7), 



(A2) 



for J/ = 23/ 1 1 . Thus our empirical model explains most of the observed variation but not all. 
We use it to estimate the step-size correction along the m,„/ 0. 1 iris trajectory. By comparison the 
plaquette slopes for the m„d ^ 0.2m.v trajectory are small. If we apply a correction according to 
a crude linear fit to these slopes, the effect on the EOS is much smaller than our statistical errors. 
For these reasons we chose to ignore the step-size error for this trajectory. 

The fit to the step size correction also allows us to estimate the error in our ability to predict the 
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FIG. 11: Plaquette slopes for the m„j ^ O-lm^ trajectory (left) and the m„j Rr: 0.2ms trajectory (right) vs. 
the chiral condensate. On the left panel the quadratic fit to the data, Eq (lAll ). is shown as well. Squares, di- 
amonds, and octagons indicate slopes determined from two different R algorithm step sizes. Fancy squares, 
fancy diamonds, and bursts indicate slopes determined by comparing R- with RHMC-algorithm measure- 
ments. 



correction. The largest error, approximately 50%, occurs at small values of the chiral condensate. 
We take this as a conservative estimate of the error in our correction throughout. 

Table HVl shows, not surprisingly, that the variation in all three gauge loop observables is corre- 
lated. With our normalization the slopes appear to be of comparable magnitude. This observation 
suggests generalizing the absolute plaquette correction to all three gauge loops. 

We also tabulate the slope of the fermion variables in Table IVl Except for the gauge contribu- 
tion, all other estimated corrections to the EOS are negligible compared with our statistical errors. 
The gauge-action correction is smaller than our statistical error in most cases; however, we have 
included it in the EOS for the niud ~ 0. 1 trajectory since, although its effect is comparable to 
the statistical error, it lowers all data points consistently. 
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p 


aniud 


anis 


Mo 


Vt^q 


Trajectory 


Vt=o 


Trajectory 


a [fm] 


*6.000 


0.0198 


0.1976 


0.8250 


1 jL 


X 4 


1800 


12^* 


500 


0.366 


*6.050 


0.0178 


0.1783 


0.8282 


123 

1 jL 


X 4 


1800 


12^ 


500 


0.334 


6.075 


0.0169 


0.1695 


0.8301 


123 

1 


X 4 


1800 






0.319 


*6.100 


0.0161 


0.1611 


0.8320 


1?3 


V 4 


1800 


12^ 


500 


0.306 


6.125 


0.0153 


0.1533 


0.8338 


1?3 


V 4 


2800 






0.293 


^^6.150 


0.0146 


0.1458 


0.8356 


123 


X 4 


3800 


12^ 


500 


0.281 


6.175 


0.0139 


0.1388 


0.8374 


1?3 

1 z. 


X 4 


3800 






0.269 


*6.200 


0.0132 


0.1322 


0.8391 


1?3 

1 z. 


X 4 

A T- 


3800 


12^ 


400 


0.258 


6.225 


0.0126 


0.126 


0.8407 


123 

1 Z. 


X 4 

A 


3800 






0.248 


^^6.250 


0.012 


0.1201 


0.8424 


123 

1 Z. 


X 4 

A 


3800 


12^ 


500 


0.238 


^^6.275 


0.0114 


0.1145 


0.8442 


1?3 
1 z. 


X 4 

A T" 


2800 


12^ 


500 


0.229 


*6.300 


0.0109 


0.1092 


0.8459 


1?3 

1 z. 


X 4 

A T" 


1800 


12^ 


2100 


0.220 


*6.350 


0.00996 


0.0996 


0.8491 


123 

1 Z. 


X 4 

A 


1800 


12^ 


2100 


0.204 


6.400 


0.00909 


0.0909 


0.8520 


1?3 

IZ. 


X 4 

A T- 


1800 






0.190 


^^6.458 


0.0082 


0.082 


0.8549 


1?3 

IZ. 


X 4 

A T- 


1800 


12^ 


2100 


0.175 


6.500 


0.00765 


0.0765 


0.8570 


1?3 

IZ. 


X 4 

A T- 


1800 






0.165 


^^6.550 


0.00705 


0.0705 


0.8593 


1?3 


X 4 

A 


1800 


20^ 


2100 


0.155 


6.600 


0.0065 


0.065 


0.8616 


123 

IZ. 


X 4 

A T" 


1800 






0.145 


*6.650 


0.00599 


0.0599 


0.8636 


123 


x4 


1800 


24^ 


2100 


0.137 


6.700 


0.00552 


0.0552 


0.8657 


123 


x4 


1800 






0.129 


*6.760 


0.005 


0.05 


0.8678 


163 


x4 


800 


243 x 64 


2100 


0.120 


*6.850 


0.00437 


0.0437 


0.8710 


163 


x4 


800 


32^ 


300 


0.109 


*7.080 


0.0031 


0.031 


0.8779 


163 


x4 


800 


403 X 96 


1500 


0.086 



TABLE I: Run parameters of the trajectory with m„^/ Rr: 0.1 at Nt = 4. The asterisk indicates paiameter 
sets for which both zero and nonzero temperature runs were performed. The columns labeled "Trajectory" 
indicate the number of thermalized trajectories. Last column shows the lattice spacing as determined from 
Eq. (E]). 
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P aniud anig uq VV^o Trajectory Vt=o Trajectory a [fm] 



*6.300 0.0109 0.1092 0.8459 12^ x 6 
*6.350 0.00996 0.0996 0.8491 12^ x 6 
6.400 0.00909 0.0909 0.8520 12^ x 6 
*6.458 0.0082 0.082 0.8549 16^ x 6 
6.500 0.00765 0.0765 0.8570 12^ x 6 
*6.550 0.00705 0.0705 0.8593 12^ x 6 
6.600 0.0065 0.065 0.8616 12^x6 
*6.650 0.00599 0.0599 0.8636 12^ x 6 
6.700 0.00552 0.0552 0.8657 12^ x 6 
*6.760 0.005 0.05 0.8678 20^ x 6 
*6.850 0.00437 0.0437 0.8710 18^ x 6 
*7.080 0.0031 0.031 0.8779 18^x6 



3100 


12^ 


2100 


0.220 


2900 


12"^ 


2100 


0.204 


2900 






0.190 


2140 


12"^ 


2100 


0.175 


2900 






0.165 


2900 


20"^ 


2100 


0.155 


2900 






0.145 


2900 


24^ 


2100 


0.137 


2900 






0.129 


1000 


24^ X 64 


2100 


0.120 


1300 


32^ 


300 


0.109 


2200 


40^ X 96 


1500 


0.086 



TABLE II: Same as Table Ubut for trajectory m„d ?a 0. 1 rris at A^, = 6. 
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p 


aniud 


arris 


"0 




Trajectory 


Vt=o 


Trajectory 


a [fm] 


*6.300 


0.0225 


0.1089 


0.8455 


12^ 


X 6 


i 3000 


12^ 


2100 


0.224 


*6.350 


0.0206 


0.1001 


0.8486 


123 

1 jL 


X 6 


i 3000 


12"^ 


2100 


0.208 


6.400 


0.01886 


0.0919 


0.8512 


123 


X 6 


i 3000 






0.193 


6.433 


0.0178 


0.087 


0.8530 


1?3 


A U 


800 






0.184 


*6.467 


0.01676 


0.0821 


0.8549 


163 


A U 


i 2200 


163 X 48 


1225 


0.176 


6.500 


0.0158 


0.0776 


0.8568 


123 


X 6 


i 3000 






0.168 


*6.525 


0.0151 


0.0744 


0.8580 


123 


X f 
A U 


i 3000 


12^ 


2100 


0.162 


6.550 


0.0145 


0.0713 


0.8592 


123 

1 Z. 


X f 
A U 


i 3000 






0.157 


*6.575 


0.0139 


0.0684 


0.8603 


123 

1 Z. 


X 6 


i 3000 


16^ 


1760 


0.152 


6.600 


0.0133 


0.0655 


0.8614 


123 

1 Z. 


X 6 


i 3000 






0.147 


-^6.650 


0.0121 


0.0602 


0.8634 


123 

1 Z. 


A U 


i 3000 


20^ 


836 


0.138 


6.700 


0.0111 


0.0553 


0.8655 


123 


x6 


i 3100 






0.130 


*6.760 


0.01 


0.05 


0.8677 


203 


x6 


i 1935 


203 X 64 


825 


0.121 


*6.850 


0.00898 


0.0439 


0.8710 


123 


x6 


i 3000 


243 


740 


0.110 


7.092 


0.00673 


0.031 


0.8781 


123 


x6 


i 3000 






0.086 


7.090 


0.0062 


0.031 


0.8782 








283 X 96 


565 





TABLE III: Same as Table Ubut for trajectory m„rf 0.2m,, at Nt = 6. The last row is a run which does not 
lie on the trajectory and was used only for zero-temperature extrapolations. 
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volume 


P 


antud 


arris 






(C) 


R 12^ X 4 


6.0 


0.0198 


0.1976 


-14(5) 


-12(6) 


-8(7) 


R 12^ X 4 


6.1 


0.0161 


0.1611 


-24(8) 


-25(9) 


-25(10) 


R 12^ X 4 


6.2 


0.0132 


0.1322 


-18(7) 


-22(8) 


-21(10) 


R 12^ X 12 


6.2 


0.0132 


0.1322 


-25(5) 


-26(5) 


-26(7) 


H 16^ X 48 


6.35 


0.00996 


0.0996 


-5(4) 


-5(6) 


-5(7) 


H 16^ X 48 


6.35 


0.0206 


0.1001 


-0.6(1.3) 


-1.0(1.6) 


-1.3(1.9) 


H 16^ X 48 


6.40 


0.00909 


0.0909 


-17(4) 


-20(5) 


-17(4) 


H 16^ X 48 


6.40 


0.01886 


0.0909 


4(3) 


5(3) 


5(4) 


R 12^ X 4 


6.458 


0.0082 


0.082 


-11(6) 


-14(8) 


-20(9) 


H 16^ X 6 


6.458 


0.0082 


0.082 


-18(4) 


-18(4) 


-13(11) 


H 16^ X 48 


6.458 


0.0082 


0.082 


-11(4) 


-11(7) 


-14(7) 


H 16^ X 6 


6.467 


0.01676 


0.0821 


-3.8(1.6) 


-5(2) 


-5(3) 


H 16^ X 48 


6.467 


0.01676 


0.0821 


2.1(1.4) 


3(2) 


4(2) 


H 16^ X 48 


6.50 


0.00765 


0.0765 


-6(5) 


-5(6) 


-4(7) 


R 12^ X 6 


6.55 


0.00705 


0.0705 


1.4(9) 


4(14) 


6(13) 


R 12^ X 12 


6.565 


0.005 


0.0484 


-36(4) 


-46(6) 


-48(6) 


R 12^ X 6 


6.65 


0.00599 


0.0599 


5(10) 


1(13) 


5(11) 


H 12^ X 4 


6.76 


0.005 


0.082 


22(11) 


19(18) 


23(16) 


R 24^ X 64 


6.76 


0.005 


0.05 


-11(2) 


-15(3) 


-10(3) 


R 20^ X 6 


6.76 


0.01 


0.05 


5(2) 


4(4) 


4(3) 


R 20^ X 64 


6.76 


0.01 


0.05 


-1.8(6) 


-2.9(4) 


-2.6(6) 


R 20^ X 64 


6.79 


0.02 


0.05 


1.24(12) 


1.12(14) 


1.23(10) 


R 20^ X 64 


6.81 


0.03 


0.05 


2.04(7) 


2.04(9) 


2.09(10) 



TABLE IV: Step-size slope do jd^ for gauge field contributions to the equation of state for a variety of 
lattice ensembles. Three operators O are tabulated. The label R indicates values determined exclusively 
from the R algorithm. The label H indicates values determined with the aid of the RHMC algorithm. 
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volume 






\ T T / J 




\ ^ duo ^ 1 


R 12^ X 4 


0.0198 


3(10) 


4(5) 


27(20) 


45(19) 


R 12^ X 4 


0.0161 


26(21) 


19(11) 


66(40) 


10(37) 


R 12^ X 4 


0.0132 


39(36) 


26(16) 


53(32) 


54(29) 


R 12^ X 12 


0.0132 


15(11) 


17(5) 


44(17) 


75(17) 


H 16'' X 48 


0.00996 


—35(14) 


-4(8) 


12(16) 


13(17) 


H 16^ X 48 


0.0206 


-2(2) 


-1(2) 


1(5) 


0(5) 


H 16^ X 48 


0.00909 


-3(9) 


3(41) 


0(18) 


0(18) 


H 16^ X 48 


0.01886 


-15(5) 


-7(4) 


7(7) 


12(7) 


R 12^ X 4 


0.0082 


1.6(1.4) 


9(6) 


27(19) 


45(19) 


H 16^ X 6 


0.0082 


5(16) 


10(8) 


36(11) 


44(12) 


H 16^ X 48 


0.0082 


-2(10) 


3(5) 






H 16^ X 6 


0.01676 


-13(4) 


-4(3) 


11(2) 


11(2) 


H 16^ X 48 


0.01676 


0.7(8) 


1.8(7) 


9(3) 


7(3) 


H 16^ X 48 


0.00765 


-5(10) 


-5(8) 


-3(15) 


15(14) 


R 12^ X 6 


0.00705 


-20(25) 


-18(23) 


-16(25) 


-6(24) 


R 12^ X 12 


0.005 


-5(9) 


18(6) 






R 12^ X 6 


0.00599 


19(12) 


33(20) 


10(22) 


5(30) 


H 12^ X 4 


0.005 


0.2(4) 


0.7(3) 


-30(28) 


-73(27) 


R 24'' X 64 


0.005 


5.4(1.8) 


1.1(1.6) 


-6.9(2.9) 


-2.7(3.3) 


R 20^ X 6 


0.01 


-1.5(1.3) 


-2.1{2.1) 


-12(5) 


-6(7) 


R 20^ X 64 


0.01 


-0.6(9) 


1.2(6) 






R 20^ X 64 


0.02 


-1.1(3) 


-1.2(2) 






R 20^ X 64 


0.03 


-0.66(11) 


-0.61(9) 







TABLE V: Step-size slope do/ dz^ for fermion contributions to the equation of state for the ensembles of 
Table JV] Four operators o are tabulated. 
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